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O 1 ABSTRACT 

04 I In the absence of any compelling physical model, cosmological systematics are of- 

ten misrepresented as statistical effects and the approach of marginalising over extra 
'r*j'' nuisance systematic parameters is used to gauge the effect of the systematic. In this 

O |. article we argue that such an approach is risky at best since the key choice of function 

I ■ can have a large effect on the resultant cosmological errors. 

As an alternative we present a functional form filling technique in which an un- 
known, residual, systematic is treated as such. Since the underlying function is un- 
known we evaluate the effect of every functional form allowed by the information 
available (either a hard boundary or some data). Using a simple toy model we in- 
troduce the formalism of functional form filling. We show that parameter errors can 
■ be dramatically affected by the choice of function in the case of marginalising over a 

systematic, but that in contrast the functional form filling approach is independent of 
the choice of basis set. 

We then apply the technique to cosmic shear shape measurement systematics and 
show that a shear calibration bias of |m(z)| < 10 _3 (1 + z) is required for a future 
all-sky photometric survey to yield unbiased cosmological parameter constraints to 
percent accuracy. 

A module associated with the work in this paper is available through the open 
source iCosmo code available at http://www.icosmo.org. 
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1 INTRODUCTION statistical power of the probes used but, most likely, by sys- 

. . . tematic effects that will be present in the data, and that are 

Cosmology is entering a formative and crucial stage, from a . , , , , , 

aj ° ° ' inherent to the methods themselves, 

mode in which data sets have been relatively small and in 

which the statistical accuracy required on parameters was The problem that we will address is how the final level of 

relatively low, into a regime in which the data sets will be systematics, at the cosmological parameter estimation stage, 

orders of magnitude larger and the statistical errors required should be treated. This problem is of relevance to all cos- 

to reveal new physics (for example modified gravity - Heav- mological probes, some examples include weak lensing and 

ens et al, 2007, Kunz & Sapone, 2007; massive neutrinos - intrinsic alignments (Heavens, Refregier & Heymans, 2000; 

Kitching et al, 2008c, Hannestad & Wong, 2007, Cooray, Crittenden et al., 2000; Brown et al., 2002; Catelan, et al., 

1999, Abazajian & Dodelson, 2003; dark energy - Albrecht 2001; Heymans & Heavens, 2003; King & Schneider, 2003; 

et al., 2006, Peacock et al., 2006) are smaller than any de- Hirata & Seljak, 2004, Bridle & King, 2007; Bridle & Ab- 

manded thus far. The ability of future experiments to con- dalla, 2007); baryon oscillations and bias (e.g. Seo & Eisen- 

strain cosmological parameters will not be limited by the stein, 2003), X-ray cluster masses and the mass-temperature 

relation (e.g. Pedersen & Dahle, 2007) to name a few. 

The general thesis we advocate in this article is that the 

* tdk@astro.ox.ac.uk standard approach to systematics, that of assuming some 
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parameterisation and fitting the extra parameters simul- 
taneously to cosmological parameters (e.g. Kitching et al., 
2008a; Bridle & King, 2007; Huterer et al., 2006; for weak 
lensing analyses), both misrepresents a systematic effect as 
a statistical signal and more importantly is not robust to 
the choice of parameterisation. 

As an alternative we will present a method, 'form filling 
functions', in which a systematic is treated as such: an un- 
known function which is present in the data. By exhaustively 
exploring the space of functions allowed by either data, sim- 
ulations or theory the effect of a systematic on cosmological 
parameter estimation - a bias in the maximum likelihood 
- can be fully characterised. This is a natural extension of 
the work presented in Amara & Refregier (2007b). We will 
present this using a simple toy model to explain and demon- 
strate the essential aspects of the formalism. We will then 
apply this technique to the problem of shape measurement 
systematics in weak lensing. 

We begin in Section 2 by categorising the different ap- 
proaches to systematics that can be taken, Section 3 in- 
troduces the parameter estimation formalism and how sys- 
tematics can be included in a number of alternative ways. 
We will then introduce a toy model that will then be used 
to introduce 'form filling functions' in Section 4 where we 
will also compare the standard approach to systematics to 
the one taken here. The application to weak lensing shape 
measurement systematics is presented in Section 5 and con- 
clusions will be discussed in Section 6. 



2 APPROACHES TO SYSTEMATIC EFFECTS 

In this Section we will discuss the problems that may be 
faced with respect to systematics and also address the pos- 
sible ways that these problems can be addressed. 

There are two scenarios in which systematic questions 
may arise. Either some data is available, from which infor- 
mation must be extracted, and the effect of systematics on 
some parameters measured must be addressed. Or one is 
planning for a future experiment and the potential impact 
of systematics on some interesting parameters must be fore- 
casted. In both cases there may be some extra data available 
that has partially measured the systematic effect, or there 
may be some hard boundary within which it is known that 
the systematic must lie - either from a theory or from simu- 
lation. When forecasting one may want to place a constraint 
on the quality of the extra data needed, or the extent of the 
hard boundary such that future measurements are robust. 
For both data fitting and forecasting there are a number of 
methods that can be employed to address the systematics 
which we review here. 

In the following we will consider a generic method for 
which the data is an observed correlation (covariance) C° hs 
which is a sum of a 'sig nal' C signal (f), which depends on a 
set of statistical (cosmological) parameters 6, and a general 
additive systematic effect C sys (that can, or cannot, depend 
on the parameters we wish to measure), so that the total 
observed signal is now 

C obs (0) = C signal (0) + C sys (0) + C noiso (0). (1) 

We have also added a benign shot noise term C n °' BC (which 



again can, or cannot, depend on the parameter(s) being mea- 
sured) . We do not claim that all systematics can be written 
this way (but most can when the data used is a correla- 
tion/covariance of quantities) - a multiplicative bias is just 
a special kind of additive term which has the same form as 
the signal but is multiplied by a systematic constant. 

We have identified three broad categories of approach 
that could taken when dealing with systematics. 

Marginalisation 

Marginalisation of systematics entails using a model, a func- 
tion containing a set of parameters a, to characterise the sys- 
tematic effect C sys — ► C sys (a) . In this case the cosmological 
parameters and the systematic parameters a are measured 
simultaneously. The extra 'nuisance' parameter errors arc 
marginalised over to arrive at the final cosmological param- 
eter errors, that now take into account the systematic. 

Marginalisation misinterprets the systematic as a sta- 
tistical signal (attempts to characterise the systematic by 
finding best fitting parameters), by reducing the estimation 
and determination of the systematic into a parameter es- 
timation problem. It would be an inappropriate statistical 
approach to estimate nuisance parameters that were known 
to have a very small degeneracy with cosmological parame- 
ters and then to claim that systematics were negligible. 

When marginalising one is immediately faced with the 
choice of model. In the absence of some underlying phys- 
ical theory one is forced to parameterise. The key choice 
of parameterisation is what makes this approach risky (at 
best); both the number of parameters and the prior (if any) 
on those parameters can dramatically affect the level of in- 
fluence that the systematic may have on cosmological pa- 
rameter estimation. One can choose either simple models, 
whose small degree of freedom may have a minimal impact 
on the cosmological parameters, but whose behaviour may 
mask the true systematic signal. Or, very flexible models; 
but one is always limited by the number of degrees of free- 
dom that can be estimated from the data, and using for 
example ^> 100 nuisance parameters to find the systematic 
error on ~ 10 cosmological parameters seems assymetric. 

In some circumstances there are physical models that 
can be called upon to model a systematic accurately, in 
this case marginalisation becomes an attractive option. 
One could also use more sophisticated techniques, such as 
Bayesian evidence, to determine which parameterisation is 
warranted given the data available. But even in such a sce- 
nario the question of whether an even more apt model is 
available, or not, would always remain and even in this case a 
residual systematic will remain (at least due to noise) which 
may contain a still unknown effect and must be treated in 
the correct way. 

Bias Formalism 

The systematic is not marginalised over in cosmological pa- 
rameter estimation, but is left in as a systematic term C sys =fc 
C sys (9). By 'leaving a systematic in' and not marginalis- 
ing over parameters the systematic is correctly identified 
as a systematic effect, albeit that the magnitude of the ef- 
fect must be correctly quantified. If a systematic is 'left in' 
then the cosmological parameter errors themselves are unaf- 
fected (in the case that the observation is not dominated by 
the systematic). The maximum likelihood value of the cos- 
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mological parameters however will always be biased by an 
amount which depends upon the true, underlying, system- 
atic signal. There have been studies of the biases that can 
be caused when a systematic is treated as such (e.g. Huterer 
& Takada, 2005; Amara & Refregier, 2007b; Kitching et al., 
2008a) although all studies have assumed some functional 
form for the underlying systematic. 

The task is then to investigate all possible functional 
forms for the systematic, that are allowed by either theory 
or data so quantifying the extent of the possible biases. In 
this case flexibility is paramount since every possible allowed 
function must be tested. This is the approach advocated in 
this article. 

By treating the systematic in this way, as a true system- 
atic effect as opposed to a statistical effect (as in marginal- 
isation), we can move away from the dilemma of choosing a 
particular parameterisation. 

Nulling 

The general approach to nulling is that the statistical signal 
used (the way in which the data is used to extract cosmo- 
logical information) can be modified in such a way that the 
systematic signal is cancelled out i.e. C obs = C signal + C sys + 
C noise - CLw = CST 1 + <?ne°w e and CZl = 0. Cosmolog- 
ical parameter estimates can then be made using this new 
statistic which by construction has minimised, or completely 
removed, the systematic effect. 

The nulling approach is a potentially powerful tool, for 
example as shown by Joachimi & Schneider (2008) this could 
be used in the removal of weak lensing intrinsic alignment 
contaminant. However when nulling, the cosmological sig- 
nal is changed in such a way that parameter constraints 
can be severely degraded. We will not address the nulling 
approach further in this article but we note that the pos- 
sibility of 'optimal' weighting (partial nulling) should exist, 
in a mean-square error sense. Nulling aims to set the bias 
due to a systematic to zero, which may be to strict since our 
true requirement is simply that the biases are sub-dominant 
to the statistical errors. 



In the next Section we will review the marginalisation 
procedure and formalise the biasing effect of leaving a sys- 
tematic in the signal. What we endorse within the context 
of the bias formalism is using the theory and data itself 
to investigate the full range of allowed functions, thus fully 
characterising the effect that a systematic may have. 



3 THE GENERAL FORMALISM 

A common approach in cosmology is to measure the signal 
of some quantity, and match this to theory in order to con- 
strain cosmological parameters. However, as we will show, 
the effect of systematics on such cosmological probes is usu- 
ally dealt with in a way which can potentially mask their 
true impact. 

Fig. 1 shows the basic situation which we will ad- 
dress. The left panel shows the observable C° bs which is 
a sum of the signal <7 sl s nal anc j some systematic plus noise 
C sys + C no,sc , equation (1). Furthermore there is some toler- 
ance envelope around the systematic (gray solid lines) which 




Figure 1. An example of the basic premise concerning the pa- 
rameter estimation methodology. The left panel shows a total 
observed correlation (solid black line), that is a sum of the sig- 
nal (dot dashed line) and systematic plus noise (dashed line). The 
systematic is known to lie within some tolerance envelope (within 
the gray solid lines). The right panel shows the observation minus 
the mean of the systematic plus noise, leaving an estimator of the 
signal (solid line) and a systematic tolerance about zero. 

represents the state of knowledge regarding the systematic. 
One can then subtract the mean C sys and C noisc from_the 
observable which results in an estimator of the signal C si s nal 

C° bS - (C SyS ) - (C noisc ) = C^gna 1 + C*y» (2) 

plus some residual systematic C sys which is centered around 
zero. In general throughout we always consider the case that 
there is some extra data that places constraints on the sys- 
tematic C sys . This is shown in the right panel of Fig. 1, the 
systematic tolerance envelope now lies about the C(x) = 
line. 

The measurement of this signal to estimate the values 
of some parameters 6 within a theory C^™ o f y (0) can be done 
in the usual way 

x\0) = E ^[c^ 1 - ctzim 2 (3) 

X 

where oc(x) is the error on the signal. Note that we will 
remove (x) (e.g crc(x) — > oc) from all equations for clarity. 
A best estimator of the parameters from the observation 6 
is defined such that d\ 2 /dO = 0. However this statistic has 
not taken into account the residual systematic effect in any 
way. 

To make predictive statements regarding parameter es- 
timation it is convenient to work with the Fisher matrix 
formalism. The Fisher matrix allows for the prediction of 
parameter errors given a specific experimental design and 
method for extracting parameters. In the case of Gaussian- 
distributed data where we assume that the error on the sig- 
nal is not a function of parameter values oc crc(0) we 
can take the covariance of the estimated values of the pa- 
rameters (Tegmark, Taylor & Heavens, 1997; Jungman et 
al., 1996; Fisher, 1935) 

cov[M,] = W - - <4»> = Fr/ (4) 

where the Fisher matrix is defined by (Tegmark, Taylor & 
Heavens, 1997; Jungman et al., 1996; Fisher, 1935) 

„ [ -2.9C OC] 

F -=2^[ a c Qe.WA- (5) 

X L J J 

The marginal errors on the parameters are given by A8i = 
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Approach 


Broadens Likelihood 


Bias Likelihood 


Problems 


Marginalise 


V 


X 


Choice of paramctcrisation 


Bias Formalism 




V 


Need to assess all allowed functions 



Table 1. A summary of the different approaches to systematic effects, showing the effect on the likelihood surface and the primary 
problem that each method encounters. ' Only in the case that the systematic is sub-dominant to the signal. 



y/ (F~ 1 )u, this is the minimum marginal error that one can 
expect for the experimental design considered (due to the 
Cramer- Rao inequality; Tegmark, Taylor & Heavens, 1997). 

3.1 Model Fitting and Marginalisation 

The marginalisation approach fits a model to the residual 
systematic and treats the systematic as an extra statistical 
effect. The model chosen for the systematic Cj^ ory (a) de- 
pends on a suite of new parameters a and on the original 
parameter set 6 where the total parameter set is given by 
<J? = {6, a). The extra parameters are then assumed to be 
part of the signal of a method. To estimate the values of the 
parameters the \ 2 statistic of equation (3) is modified to 



Xtotai(0,o) = 



JO 



E< 

X 

E 



gnal 



^signal /n\ 
° theory \") 



c sys 

theory 



(a)f + 



-2 



/~fsys 
theory 



(a)] 2 



(6) 



where the total \ 2 is minimised to find the best estimator 
of the parameters. The likelihood functions for the cosmo- 
logical parameters are found by marginalising the combined 
likelihood p(6, a) over the new parameters 



P(0) 



dap{6, a) 



(7) 



The new Fisher matrix for the total parameter set <]? be- 
comes a combination of the cosmological Fisher matrix F ee , 
the derivatives of the likelihood with respect to the cosmo- 
logical parameters and the systematic parameters F ea and 
the systematic parameters with themselves F aa 

! F ee F e 



F c 



F c 



(8) 



Where the individual terms are given by 
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{^signal " 



fl/^signal o^signal 
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L "~'theory "'-'theory 
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(^signal " 



"'-'theory "'-'theory 



da i 
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— 2 theory theory 



(9) 



here we have assumed that the errors are uncorrelated and 
do not depend on the parameters. 

The predicted cosmological parameter errors now in- 
cluding the effect of the systematic are given by A6i = 



( see Appendix A for a more detailed expres- 
sion). The cosmological parameter errors are increased due 
to the degeneracy between the cosmological and systematic 
parameters (included by the F 6a terms). The tolerance en- 
velope around the residual systematic (Fig. 1) acts as a prior 
on the systematic parameters in the chosen model. 

3.2 The Bias Formalism 

The bias formalism treats a systematic as such, by not sta- 
tistically marginalising over any extra parameters within a 
model. Instead the systematic simply adds an extra system- 
atic function to the signal. By doing this a bias is introduced 
in the maximum likelihood value of the parameters with re- 
spect to the true underlying values 



b[Oi] = (9i) - {Of™}. 



(10) 



When marginalising the choice lies in the suite of parame- 
ters, and the function, chosen. Here there is a similar choice, 
one must assume that the systematic has some functional 
form Cf^ ction . To estimate the values of the cosmological 
parameters 8, the \ 2 statistic of equation (3) is modified to 
include the assumed systematic 



x 2 (0) = ]rac 2 [c^ + a 



sys 

function 



^signal /a\]2 
U theory I. a )\ ■ 



(11) 



The estimate of the parameter values will now be biased 
but the marginal error on the parameters will remain the 
same (the caveat here that this is only the case when the 
systematic is smaller than the signal). 

It can be shown (Taylor et al., 2007; Amara & Refregier, 
2007b; Kim et al., 2004) that, with the assumption of Gaus- 
sian likelihoods, the predicted bias in a parameter due to an 
uncorrected systematic is given by 

b[0i] = (F~ 1 )ijBj 



(12) 



where 



p. . _ 



= E 

X 

= E 



o^signal o^signal 
-2 "^theory "°theory 



dOi 



a C ^function" 



£\s~isignal 
theory 



d0< 



(13) 



To recap Sections 2 and 3, Table 1 summarises the effect 
on the likelihood surface and the primary problem encoun- 
tered by each systematic approach. 

In cases where the systematic affects the error on the 
signal o c = crc(C sys ) (which is almost always) then leaving 
a systematic in the signal can cause a bias and increase 
the marginal errors. However the increase in marginal errors 
is negligible for systematics that have an amplitude which 
is much less than the signal, and cause biases that are < 
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10cr (Amara & Refregier, 2007). Equation (12) is also an 
approximation for the case of small biases, if the bias is 
large relative to the marginal error then the curvature of 
the likelihood surface will have varied substantially from the 
Gaussian approximation. In such cases one could go to a 
higher order in the Taylor expansion used to derive equation 
(12), or calculate the full likelihood. 



3.3 A Simple Example 

To review the general formalism described thus far, and for 
use in subsequent Sections, we will here introduce a simple 
model. The toy model we will consider is shown in the right 
hand panel of Fig. 1. Referring to equation (1) the signal is 
given by a simple polynomial expansion 



signal 



; (x) = aa + aix + (-0.45)x 2 + (0.05)x 3 



(14) 



J example \ 

where the statistical parameters we are concerned with (with 
fiducial, true, values) are the parameters ao — 1.0 and 
a± = 1.0. We assume that the observed signal is measured 
with an error of ac(x) = 0.25. We include the x 2 and x 3 
terms so that the problem is slighly more realistic, in that 
there is an extra behaviour in the signal that we do not 
wish to constrain but may effect the determination of the 
parameters of interest. 

Fig. 2 shows the simple model signal with some Gaus- 
sian distributed data points. We then calculate the two- 
parameter marginal errors using equation (3) where we use 
the model as in equation (14). Table 2 shows the measured 
marginal errors from this simple mock data, and compares 
these with the expected marginal errors calculated using the 
Fisher matrix (equation 5). 

We then introduce a simple systematic into the model 
by assuming a simple function (note this does not have a 
mean of zero, but could fit into some boundary centered on 
C(x) = 0) 



example 



(a;) = -0.2 + 0.15a;- 0.01a; 2 



(15) 



By recalculating the likelihood and including this systematic 
function, using equation (11), we find that the most likely 
value of oo and ai is biased and yet the marginal errors 
remain the same. We compare this bias to the prediction 
made using equation (12) in Table 2. 

To test the method of marginalising over data we now 
reset the systematic (throw away equation, 15) and intro- 
duce a simple systematic which is measured using some 
mock data. Fig. 3 shows the model signal with some ex- 
tra systematic data with a mean of zero and a scatter of 
(Tcsr- =0.5. We then introduce a simple systematic model 
parameterised by a new parameter so 



'-'theory example \ X ) ~ S ° 



(16) 



and fit the model (equation 14) and the systematic to the 
data simultaneously, as described in equation (6) . The total 
likelihood p(ao,ai, so) is then marginalised over so. Fig. 3 
shows that when this extra systematic is marginalised over 
the constraint on ao is affected the most since so has the 
same functional form as this parameter, and so there exists 
a large degeneracy between ao and so . The exact degeneracy 
is broken by the extra data available for the systematic. In 
Table 2 we show the increase in the marginal error on the 





T).B 0.9 



8 10 



Figure 2. The left panel shows the toy model signal, the fiducial 
model is shown in red (dark gray) and we have added some sim- 
ple Gaussian distributed data points about this central model. 
We also show an example systematic in (light) gray defined in 
equation (15). The right panel shows the two-parameter 1-tr error 
contours without a systematic (black) and including the example 
systematic (red/dark gray). The fiducial model is marked by +. 



No Systematics 






Parameter 


Measured Error 


Expected Error 


ao 


0.100 


0.099 


ai 


0.011 


0.011 


Bias Method 


Parameter 


Measured Bias 


Expected Bias 


ao 


0.003 


0.002 


«i 


0.040 


0.040 


Marginalising Method 


Parameter 


Measured Error 


Expected Error 


ao 


0.120 


0.119 


<ii 


0.012 


0.013 



Table 2. This table compares the Fisher matrix predictions with 
values found using some simple mock data described in Section 
3.3. This is not a comparison of the systematic methods them- 
selves, which will be done in Section 4.3. The upper table shows 
the marginal errors found using the data shown in Fig. 2 for the 
parameters ao and a\ defined in equation (14). These errors are 
compared to what is expected from the Fisher matrix, equation 
(5). The middle table shows the measured bias when a simple 
systematic is added (equation 15) and compares this to the ex- 
pected bias calculated using equation (12). The lower table shows 
the increased marginalised errors on ao and ai when a simple pa- 
rameterised systematic model is marginalised over. 



parameters ao and a\ when we marginalise over the extra 
systematic parameter so- 

We have now introduced the basic formalism and shown 
that this can be applied to a simple example that yields 
results which are in good agreement with the Fisher matrix 
predictions. 



4 FORM FILLING FUNCTIONS 

For the remainder of this article we will present an alter- 
native to marginalisation by advocating the bias formalism 
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Figure 3. The left panel shows the model signal, the fiducial 
model is shown in red (dark gray) and we have added some 
Gaussian distributed data points about this central model. We 
also show an example data-driven systematic blue (light gray) 
points about C(x) = with a variance of crpsys = 0.5. The 
right panel shows the two-parameter 1-cr error contours with 
(green/light gray) and without (black) marginalising over the sys- 
tematic model. The red lines show the marginal error with no 
extra systematic data, a complete degeneracy between arj ans so- 
The fiducial model is marked by +. 
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Figure 4. Representing the two possible systematic constraints, 
either from theory or from data. A hard boundary (left panel solid 
gray lines) may be defined within which the systematic must lie, 
or some data may provide a measurement of the level of system- 
atic (right panel blue/light gray error bars). The black data points 
are mock data with a Gaussian distribution using the model de- 
scribed in Section 3.3, the red (dark gray) line is the fiducial 
model. 



for dealing with systematics, outlined in Section 3.2. The 
issue with which one is now faced is what function to choose 
for the residual systematic. To investigate the full extent of 
possible biases, allowed by the tolerance on the systematic, 
all allowed functions must be addressed in some way. 

In Fig. 4 we use an extension of the simple model, out- 
lined in Section 3.3, to introduce the concept of two different 
forms of systematic tolerance. The tolerance envelope could 
have a hard boundary (e.g. defined by a theory which states 
that "the systematic must lie within this boundary"). Or 
the systematic could be defined by some extra data that has 
partially measured the magnitude of the systematic. We will 
refer to these two scenarios as the "hard bound" and "data 
bound" respectively. 

To fully assess the level of bias every functional form 
allowed by the systematic tolerance envelope needs to be 



tested. For the hard boundary we want to find every func- 
tion that can be drawn within the hard boundary. For the 
data bound every function can be weighted with respect to 
the data itself. Here we introduce the concept of 'form filling 
functions' which are a set of functions that should exhaus- 
tively fill the space of possible functions allowed by some 
tolerance envelope. 

Consider the hard bound in Fig. 4, in the bias approach 
we want to find every function that will fit within this toler- 
ance envelope. To do this we consider functions in the most 
general form as expansions in some arbitrary basis set 

JV 

f(x; {a n }, {b„}) = 2J a n tpn{x) + b n cj> n (x) (17) 

n=l 

where a n , b n 6 5ft and xj) n {x) and 4> n (x) form some arbi- 
trary basis functions. To choose the basis set we impose the 
following conditions 

• The basis set must be complete in the range of x we are 
considering i.e. all functions f(x) must be expressible as an 
expansion in the basis. 

• The basis functions must be orthogonal. 

• The functions must be boundable i.e. the basis set must 
be able to be manipulated such that every function described 
within some bounded region (tolerance envelope) can be 
drawn. 

The first condition is necessary. The second condition is 
make some calculations more straightforward though is not 
strictly necessary. The third condition is merely desirable - 
one can imagine having a basis set in which non-bounded 
functions are allowed, but when using such expansions one 
would have remove these stray functions. As an aside we 
note that over-completeness is not a problem as long as the 
basis set is in fact complete (we do not care if we sample 
a function multiple times, as long as we sample it at least 
once)^. 

Algorithm. Equation (17) represents every function, how- 
ever we are only concerned here with defining all functions 
within some bounded region. To begin we will show how an 
interval \a n \ < Q can be defined such that every function 
between |/(a;)| < 1 can be drawn. For an orthonormal basis 
set the coefficients needed to draw a function f(x) can be 
expressed as 

/ dxw(x)f(x)(/> m (x) (18) 

JR 

where R is some interval over which the basis set is complete 
and w(x) is a weight function upon which the basis set is 
complete the constant A (A mn is a diagonal) is calulated in 
general using 

A^m = / dxw(x)(p n (x)4> m (x) (19) 

J R 

for all the basis sets we consider A nm = AS„ m - Now we can 
write a general expression that provides a limit on a n . Using 

t Orthogonal basis sets are never over-complete so the second 
condition means we shouldn't be in this situation for our func- 
tional form filling algorithm. 
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Figure 5. An example of a bounded area (thick black lines), and a random sampling of 15 functions shown by thin coloured (gray and 
black) lines. Every function is defined, using equation (22) so that it must fit within the bounded area. In the left panel we used the 
Chebyshev basis set, in the central panel we use the Fourier basis set and in the right panel we use the tophat basis set (binning). For 
all basis sets the maximum order is N = 15. The functions are defined by uniformly and randomly sampling the coefficient space {ai, 
015}. This Figure is ment simply as a example of the type of function that can be drawn not as a measure of wether the method 
succeeds in drawing all functions. The success of the method in drawing all functions is shown in an extensive fashion in Appendix B. 



the triangle inequality we can write 

\a n \ <a[ dx\w(x)\\<j> n (x)\\f(x)\ 

J R 

and given that |/(a;)| < 1 we have an expression for \a 
\a„\ < A 



(20) 



/ dx\w(x)\\cl> n (x)\\f(x)\<A[ dx\w(x)\\4>„( x )\ = 

J R J R 

(21) 

So using equation (17) and limiting the coefficient values 
to \a n \ < Q (and similarly for b n ) every function with the 
region |/(a;)| < 1 can be drawnt. 

We now define an arbitrary 'bound function' B(x) 
which describes a hard boundary (in Fig. 4 for example) 
where at any given point in x the systematic functional form 
must lie in the region —B{x) < f(x) < B(x). Equation (17) 
is simply modified to include this arbitrary boundary 

JV 

f(x; {a n }, {b n }) = B(x) a n i> n {x) + b n <j>„(x). (22) 

n 

Now, if the basis set is complete, N = oo and \a n \ < Q 
(similar for b n ) equation (22) represents every possible func- 
tion that can be drawn within the hard bound - and each 
function could yield a different bias. Note however that this 
statement says nothing about the probability that a partic- 
ular function will be drawn. 

An important caveat to this is that the total set of func- 
tions drawn using this algorithm is not bounded only some 
subset of the functions is. However a 'clean' subset of func- 
tions with |/(a;)| < B(x) can easily be drawn by removing 
any function for which > B{x) at any x. 

In practice where N < oo the task of drawing all possi- 
ble functions becomes a numerical/computational problem. 
For a given basis set the fundamental quantities that de- 
scribe each and every function are the coefficients a n and 

■f Equation (21) is strictly only true for Riemann integrable func- 
tions, however essentially all bounded functions satisfy this con- 
straint - in particular all class C° (smooth) functions, and all 
step functions. Examples of the type of (very peculiar) function 
that are not Riemann integrable are Dirichlet's function and the 
Smith- Volterra-Cantor set. 



b n . The task then is to explore the coefficient parameter 
space {\a n \ < Q} in an exhaustive a manner as possible. 

As a first attempt the approach taken in this article is 
to randomly and uniformly sample the coefficient space {ai , 
ajv} for \a n \ < Q. The free parameters in this approach, 
given a basis set, are the maximum order investigated JV 
and the number of random samples in the space {a n } that 
are chosen. We leave a more sophisticated Monte-Carlo for- 
mulation of this problem for future work. 

When truncating the series we will be missing some 
highly oscillatory functions (for the basis sets considered) 
but we show in Appendix B that one can always make a 
definitive statement about the fraction of all functional be- 
haviour sampled down to some scale. These free-form func- 
tions are regularised by the truncation of the basis set and 
by the bound function. Throughout the remainder of this 
Section we use the numerical order and function number in- 
vestigated in Appendix B. 

Basis Sets. Throughout we will principally consider three 
different basis sets. These are 

• Chebyshev polynomials T n (x). These functions form a 
complete basis set for — 1 < x < 1, and are bounded by the 
region |/(a;)| < 1. To map these functions onto an arbitrary 
arrange a variable transformation can be applied such that 



ip n (x) = cos(n arccos(a;)) = T n 



2x - 3^ mi ii Xjj 



, (23) 



and 4> n {x) = for all n. 

• Fourier series. The Fourier series is a complete basis set 
in the range — n < x < n and the functions are bounded in 
the region |/(x)j < 1 for the basis set 



1pn(x) 
<f>n(x) ■ 



2 
1 



irx 



nil [i 



3-max 

nx — 



X 



mm 



(24) 



The cosine or sine of a real number are |cos(a;)| < 1 and 
sin(x-)| < 1. 

• Tophat functions (binning). This uses a tophat func- 
tional form and is meant to be analogous to binning the 



© 0000 RAS, MNRAS 000, 000-000 



8 T. D. Kitchmg et al. 



rr-range. The maximum order in equation (17) TV here refers 
to the number of bins where for the n th bin i>(x) is either 
zero or one depending on whether x lies within the bin 

, , , fl V (x„ - Ax/2) < x < (x„ + Ax/2) 

ll>n{x) = < Q 

= H(x-^)-H(x+^) (25) 

where x n = x miri + (n—l)Ax and Ax = (:r max — x min )/(N— 1) 
is the bin width. 4>n(x) = for all n. 

Table 3 summarises some of the basis set properties includ- 
ing the coefficient intervals over which a complete (sub)set 
of functions with \f(x)\ < 1 can be drawn. 

One may be concerned that there will always be some 
x at which a boundary-touching function cannot be drawn. 
For every finite maximum index N in the sum of equation 
(22), one can find a value of x in R such that f(x) cannot be 
1 at a;. The key realisation that we stress here is that this is 
a actually a statement about the scale upon which the space 
of functions is complete. The larger the maximum order N, 
the smaller one can find an e with xp n (x±e) — 1 for some n, 
i.e. the concern is reduced to an issue of resolution because 
if e <C s then every function complete down to some scale s 
can be drawn. We show this in Appendix B. 

4.1 Hard Bound 

Fig. 5 shows an example of a hard boundary and a random 
assortment of functions (a random, uniform, sampling of 
{oi, ajv}; and {&i, 6jv} for the Fourier basis set) for 
the Chebyshev, Fourier and tophat basis sets. As the order 
and number of functions is increased the functional forms 
begin to completely fill in the bounded area (hence "form 
filling functions"). 

It can be seen even at this stage that the tophat basis 
set is not efficient at filling the bounded area. This is investi- 
gated further in Appendix B where we show that whilst all of 
the basis sets considered can fill any desired bounded region 
the Chebyshev and Fourier basis are many orders of magni- 
tude more efficient in terms of computational time than the 
tophat functions (binning); we discuss computational time 
in Appendix C. This is a result of the restrictive step-like na- 
ture of the functions that binning imposes requiring a high 
order (number of bins) to characterise particular (smooth) 
functional behaviour. In addition to being inefficient, the 
tophat basis is not differentiable at the bin boundaries and 
as such could be construed as being un-physical, although in 
some special circumstances (e.g. photometric redshifts where 
a filter may have a tophat function in wavelength) a tophat 
basis set may be needed. 

The key feature of the hard bound is that within the 
boundary all functions are given equal weight i.e. the prob- 
ability that any given function is the 'true' systematic func- 
tional form is the same for all functions. Out of the full range 
of possible biases given a hard boundary, there should exist 
a maximum bias - because the space of functions and hence 
the range of biases is limited. We show that this is the case 
in Appendix D. The quantity of interest is therefore this 
maximum bias that is allowed by the functions that can 
be drawn within the hard boundary. In the data boundary 
case, Section 4.2, there is a maximum bias for each subset 



of functions that give the same weight (with respect to the 
data) . 

One intuitively expects that the systematic function 
that introduces the largest bias should be the one that most 
closely matches the signal term containing each parameter. 
In the case where there are degeneracies between parame- 
ters the worst function is some combination of the signals 
sensitivity to the parameters. Here we will use the simple 
example hard bound from Section 3.3 to demonstrate that 
in this case a maximum bias exists, and that this maximum 
is stable with respect to basis set. 

Fig. 6 shows the maximum bias on ao as a function of 
the number of random realisations of the coefficient space for 
the Chebyshev, Fourier and tophat systematic basis sets, for 
all sets we consider a maximum order of 100. It can be seen 
that as soon as the "worst function" is found the maximum 
bias becomes constant and stable for the Chebyshev and 
Fourier basis sets. In this case the worst function is simply 
C sys (a;) = x since this is the function that affects the signal 
(C sls (x) — ao + a\x + constant) the most through the ef- 
fect on a\ (the worst function is a linear combination of the 
response of each individual parameter in the marginalised 
case, see Appendix D). There is an exact degeneracy since 
a function of the form —f(x) will cause a bias of equal mag- 
nitude but opposite sign to f(x), so C sys (:r) = —x in this 
case will also cause the same absolute bias (see Appendix 
D). For the top hat basis set it is very unlikely to find 
this particular worst function so the bias does not find the 
maximum even after 100 realisations. The Chebyshev ba- 
sis set finds this function after only a few realisations since 
rpi (x) = Ti (x) — x is one of the basis functions of the Cheby- 
shev expansion. 

Fig. 6 is meant simply as an example of the type of con- 
vergence test that could be performed given a realistic ap- 
plication. In Appendix B we outline a diagnostic mechanism 
that can gauge whether a bounded area has been completely 
filled with functions. 



4.2 Data Bound 

If the systematic has been partially measured using data 
then there is no hard boundary within which all functions 
must lie, and every function is still allowed to some degree. 
The issue which must be addressed in this case is how each 
function should be weighted given the data available - a 
function is more likely to be the 'true' systematic signal if 
it is a good fit to the systematic data. 

For the case of a data bound we propose a similar ap- 
proach to the hard bound. We want to explore all possible 
functions, to do this we express an arbitrary function using 
an expansion as in equation (22). By choosing a complete 
basis set (e.g. Chebyshev, Fourier) the full space of func- 
tions can be explored by exhaustively exploring the space of 
coefficients, as described in Section 4.1. 

In the following we will concatenate the two basis sets 
{a n } and {&„} for clarity, however note that for the Fourier 
basis set two sets of coefficients are needed. 

The critical difference in this case is that for each 
function f(x; {a n }) that is drawn we can assign a weight 
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Basis Set 


Basis Functions 


Orthogonal Weight w(x) 


Orthogonal Constant A 


Interval R 


Coefficient Interval Q 


Chebyshev 

Fourier 

Tophat 


T n (x) = cos(ra arccos(rr)) 
i for n = 0, cos(?i:r) & sin(nx) 


(1 - x 2 )~i 
1 
1 


1 for n = 0; ^ for n + 

7T 7T ' 
77 

1 


[-1-1] 

[ — 7T, 7r] 

[-1,1] 


1 for n = 0; A for n ^ 

2 for n = 0; * for n ± 

I 



Table 3. This table lists some constants and functions associated with the basis sets used in this article, Chebyshev, Fourier and tophat 
functions. Some of these constants are defined before and used in equation (21). The coefficient interval is such that for \a n \ < Q all 
function in with < 1 can be drawn. 




Figure 6. Using the simple example from Section 3.3. The max- 
imum bias on arj as a function of the number of realisations of 
the systematic basis set's coefficients. Red (dark gray, upper line) 
shows the maximum bias found using the Chebyshev basis set, 
green (light, middle line) gray shows the maximum bias for the 
Fourier basis set and blue (darkest gray, lower line) for the tophat 
basis set. 



iy({a„}) using the x 2 statistic 



E 



[/(a; {a n }) - C^f 



(26) 



This quantifies how good a fit the function (the values of 
a n and b n , and some basis set) is to the residual systematic 
data. 

However since we aim to exhaustively try every function 
(complete down to some scale) there will exist some func- 
tions that fit exactly through the data points. This causes 
a problem when using the x 2 statistic as a weight for such 
functions will have a x 2 = 0, and hence be given a probabil- 
ity P = 1 that such a function is likely, but such a conclusion 
has not taken account of the error bars on the systematic 
data. 

In this simple example we have subtracted the mean 
of the systematic signal so data should be scattered about 
C(x) = 0. In the case of Gaussian distributed data the scat- 
ter of the data points should be proportional to the error 
bar on each data point. If the observation could be repeated 
then the data points would be scattered in the same statisti- 
cal manner about C(x) = but have different actual values. 



This is very similar to the familiar sample variance; in cos- 
mology we are used to a special kind of sample variance that 
we call cosmic variance in which the likelihood of the data, 
given only one realisation of our Universe, must be taken 
into account. 

To take into account this sample variance effect we must 
consider the likelihood of the residual data. In Appendix E 
we show how for Gaussian distributed data the probability 
of a function f(x; {a„}), given some set of observations, can 
be written as a sum over x, where at each point there is 
some data with an error bar acy as 



mp[/({a n })] oc 



E 



/(z;K}) 2 



ln(<7c=y= V^7r) 



(27) 



where we have translated the notation of equation (91) to 
reflect that of equation (26). For any given set of data the 
second term in equation (27) is a benign additive constant 
so the weight that we assign each function is 



W({a n }) = J2 



f(x;{a n }) 2 



(28) 



We stress that this is for Gaussian data with a known mean 
of zero such a function is uniquely described by the variance. 

In general data may not be exactly centered at zero 
or be Gaussian distributed. As an alternative to the ana- 
lytic procedure of marginalising over the data one could cre- 
ate Monte-Carlo realisations. For each realisation one could 
measure the \ 2 °f the fit of a given function to the data 
from equation (26) and assign a weight as the average over 
all realisations 



^(K}) = i(E- ;K}) 



(29) 



realisations 



where the extra factor of 1/2 converts the average x 2 to a 
log- likelihood similar to equation (28). 

The bound function in equation (22) represents a hard 
prior in this case. One should choose a bound function is 
much larger than the scatter in the systematic data B(x) 3> 
a B y B {x). Any functions that deviate from the data by a large 
amount will be down-weighted by the poor fit even though 
they may yield a large bias, so as long as the boundary 
within which functions are considered is B(x) > 3<j aya (x) 
away from the data then any results should be robust. 

For each function f(x; {a n }) one now has an associated 
bias and a weight. Now consider a particular bias in some 
parameter: there exists a set of systematic functions that 
could yield this same bias and from that set there must 
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Figure 7. The three scatter plots show the bias in the parameter 
ao caused by fitting functions through the data bounded system- 
atic shown in Fig 4 against the weight given to the function with 
respect to the systematic data points given in equation (28). Each 
point represents a function. The colours correspond to the differ- 
ent basis sets considered blue (darkest gray, bottom left) is tophat 
functions, red (lighter gray, top left) is Chebyshev functions and 
green (lightest gray, top right) is Fourier functions. Also shown 
is the likelihood of the bias in ao for each basis set, found by 
measuring the minimum weight for a particular bias - the low- 
est extent of the scattered points in the other plots — and using 
equation (30). 

exist a function which is the best fit to the data. If each 
best-fitting function for every bias can be found then one is 
left with a robust weight for each bias and a measure of the 
relative probability allowed by the data 

p(bi) oc exp(- min[W(6 i; {a„})]). (30) 

min[W(bi; {a n })] is the minimum weight (equation 28) from 
the space of functions defined by the coefficients {a n } that 
yields the bias hi. 

Another way of putting this is that for a given weight 
there exists a maximum and minimum bias. We show in 
Appendix D that, within the Fisher matrix approximation 
of equation (12), that there does indeed exist a maximum 
and minimum bias for each weight. 

Fig. 7 shows the weight (equation 28) for each func- 
tion drawn from the Chebyshev, Fourier and tophat basis 
sets against the bias induced by the function using the data 
bound toy model of Fig. 4. It can be seen from this Figure 
that for any given bias there exists a minimum weight that 
can be achieved by the functions giving that bias. In the 
right bottom panel of Fig. 7 we have found the minimum 
weight for each bias and converted this into a likelihood us- 
ing equation (30) showing that this procedure is robust to 
the choice of basis set used (though the tophat basis set is 
far less efficient than Chebyshev and Fourier, see Appendix 
B for more details). 

One can of course extend the formalism introduced here 
to an arbitrary number of dimensions. To extend this to two 
dimensions we have modified the simple example to have a 
much smaller data bound, since (as can be seen in Fig. 7) the 
fiducial scatter in the toy model systematic data introduces 




x a Q 

Figure 8. The left hand panel shows the modified simple ex- 
ample, see Section 3.3, where we have made a more tightly con- 
strained Gaussian residual systematic with crosys = 0.05 for clar- 
ity, this can be compared with the right hand panel in Fig. 4. The 
right hand panel shows the 1 — a statistical constraints on ao and 
a\ (black solid line). In addition we show the systematic proba- 
bility contours for the residual data bound in the left hand panel. 
Shown are the 1 — <r systematic contours using both a Chebyshev 
(red, dark gray line) and a Fourier (green, light gray line) basis 
set for functional form filling. 

large biases in the measured parameters. Fig. 8 shows the 
joint 1 — a statistical constraint on ao and ai for the mock 
data shown. On this same plot we show the systematic 1 — a 
contours from the residual systematic - within the system- 
atic contour there is a > 68% probability that the statistical 
maximum likelihood is biased. 

We again show results for two different basis sets and 
show that the systematic contours are again independent 
of the basis set chosen. We have not used the tophat basis 
because it is not an efficient basis set (in terms of computa- 
tional time, see Appendix B) to use even for a one dimen- 
sional analysis. 

In a real application one would hope to show such plots 
with the systematic contours lying within the statistical 
ones, but for illustrative purposes we have shown a dominant 
systematic here. A hard bound in this case would represent 
a tophat contour in probability where, within the contour, 
the probability equals unity and outside the contour the 
probability is zero. 

In the marginalisation approach one could also draw 
two sets of contours, but both would be statistical: one that 
has no systematic and the other in which a systematic has 
been included. Fig. 8 represents one of the key recommen- 
dations of this article, that in future we must not only show 
statistical contours for cosmological parameters but also sys- 
tematic probability contours. Here we have shown one way 
of obtaining a robust estimate of such systematic probability 
contours. 



4.3 Marginalisation vs. Bias 

The bias and the marginal error are inter-related via the 
mean square error (MSE) which is defined, for a parameter 

ai, as 

MSE = <j 2 (ai) + b 2 (a,i) (31) 
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were marginalising a choice of parameterisation (basis and 
order) would have to made in which case the MSE conclu- 
sion would be dependent on this choice. 
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Figure 9. The weight (fit to the systematic data) against the bias 
in the parameter ao for the toy data bounded systematic shown in 
Fig 4. We trucate the Chebyshev and Fourier basis sets at N = 1 
(black dots) TV = 10 (red, green — light gray dots) and TV = 50 
(blue, dark gray dots). For each order we make 500 realisations 
of the basis set. 



where a(aj) is the marginal error on en and b(a,i) is the 
bias. So both marginalisation over systematic parameters 
and functional form filling increase the MSE through the 
marginal error and bias respectively. 

In Appendix A we show that the marginalising ap- 
proach and the functional form filling approach are different 
however there is a subtelty. If the parameterisation and or- 
der are the same (i.e. a truncated basis set is used) then the 
MSE recovered from marginalising or considering the bias is 
the same. The conclusions from Appendix A are 

• If the functional form of the systematic is known then 
the degradation in the MSE as a result of marginalising or 
functional form filling (bias) is the same. 

• If the functional form of the systematic is unknown then 
the MSE from marginalising will tend to underestimate the 
true systematic error and is in general not equal to the MSE 
from functional form filling. 

• Given that the marginalising necessarily truncates the 
basis set there are always some functions that marginalising 
cannot assess. 

Given that the MSE, given a particular functional sub- 
set is the same for bias and marginalisation the key differ- 
ence between marginalisation and functional form filling is 
that in marginalisation case the functional space is trun- 
cated by the choice of parameterisation and ultimately by 
the number of data points. 

Looking back at Fig. 7 the functional form filling tech- 
nique fills out the Weight-MSE (bias) bounded region by 
sampling every function down to some scale. To illustrate the 
way in which marginalisation cannot sample the full space 
of functions we have reproduced these plots but using vari- 
ous truncated basis sets. Note that the MSE from the bias 
and marginalisation is the same given a parameterisation as 
shown in Appendix A - marginalising is like functional form 
filling but with a very restricted class of function. 

Fig. 9 shows the weight (fit to data) vs. bias using the 
simple toy model for the Chebyshev and Fourier truncated 
basis sets for 500 realisations of each truncated set. We trun- 
cate the expansion at N = 1, 10 and 50. It can be seen that if 
TV = 1 is used the space of functions is very restricted, and 
as the order is increased the space of functions increases. 
These results are also in resonance with Appedix B. If one 



Discussion 

We take this opportunity to discuss the difference be- 
tween marginalising and functional form filling. Functional 
form filling is not the same as marginalising over a very large 
parameter set on a qualtitative level: we are not finding the 
best fit values of the parameters but rather using each pa- 
rameter combination simply as a prescription that yields a 
particular bias. As the freedom in the functions increases 
(more parameters) the functional form filling approaches a 
stable regime in the results it gives. Functional form filling 
yields the probability that the maximum likelihood value 
is biased by some amount, whereas marginalisation yields 
a probability that the cosmological parameters take some 
value jointly with some values of the nuisance parameters. 

One could marginalise over a very large number of pa- 
rameters but in this case the joint covariance matrix will 
at some point necessarily become computationally singular 
as parameters are introduced can cannot be constrained by 
the data. One way to consistently include more parameters 
than data points TV is to add a prior such that any extra pa- 
rameters were constrained P{A\D) = J dBP(A, B\D)P{B) 
where A — TV ansd (A + B) > TV; but this requires adding 
the prior on P(B) which requires justification. 

In the lack of any external prior the number of pa- 
rameters that can be used in marginalising is necessarily 
truncated at a low order since once the number of free pa- 
rameters becomes larger than the number of data points 
(for the data bound) the parameter fitting methodology be- 
comes ill-posed (e.g. Sivia, 1996). In contrast the functional 
form filling technique could use an infinite basis set expan- 
sion (only computational, and physical, resources prohibit 
this) since we use the basis set simply as a prescription for 
drawing functions through the systematic - one could even 
draw these by hand if you were sure that you could draw a 
complete set of functions. 

One may be concerned that information is being ignored 
or disgarded. Quite the opposite from discarding information 
present in data we assign a weight to every possible function 
with respect to the data, in this sense we throw nothing 
away - every function has a weight and a bias. In constrast 
when marginalising the functional space assessed is limited 
by the number of data points and as such a very large class 
of functions are never considered - if a basis set is truncated 
these are usually highly oscillatory functions. 

The bias formalism becomes preferable if there is un- 
certainty over what the functional form is (which is most, if 
not all of the time). 

In Appendix A we show that for the same basis set and 
prior functional form filling and marginalisation produce the 
same results in a mean-square-error sense. Hence in the limit 
of marginalising over functions, with suitable priors, the two 
approaches should produce the same results - in this article 
we focus on bias functional form filling. 
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4.4 Summary 

We have presented an alogorithm by which any function 
within a bounded region can be drawn this is summarised 

as 

• Chose a complete orthogonal basis set. We recommend 
Chebyshev polynomials because of their ease of calculation 
and computational efficiency (see Appendix B). 

• Use equation (21) to calculate the interval a n € [-Q, Q] 
over which the coefficient must be sampled. 

• Define the bounded region B(x) within which functions 
must be drawn. 

• Define a scale Ax upon which the functional space must 
be complete. 

• Randomly and uniformly sample from the space of co- 
efficients using a maximum order N and number of reali- 
sations such that the functional space is fully filled - the 
diagnostic tools from Appendix B can be used to guage the 
level of completness. 

By defining the hard and data bounds we have now pre- 
sented all the tools needed to correctly assess a systematic 
given some prior knowledge of the magnitude of the effect, 
either an external data set or some theoretical knowledge. In 
the hard boundary case every function is given equal weight 
and as such a maximum bias should exist. In the data bound- 
ary case a probability can be assigned to each bias. 

We emphasise here that this method requires there to 
be at least some information at every data point associated 
with the signal - either a hard boundary or some systematic 
data. If there were no constraint the ranges of biases could 
become unbounded. 

Our proposal is that when measuring some cosmologi- 
cal parameters these techniques can be used to augment the 
statistical marginal error contours: some cosmological pa- 
rameters are measured and about their maximum likelihood 
point are drawn some marginalised statistical error contours 
and in addition: 

• If a theory or simulation provides a hard boundary to 
some systematic then the maximum bias will define a sys- 
tematic contour that can be drawn in within which the max- 
imum likelihood could be biased. 

• If some data is provided that measures the systematic 
then a further set of systematic contours can be drawn which 
will show the probability that the maximum likelihood point 
is biased by any particular amount. 

The goal for any experiment is to control systematic effects 
to such a degree that any systematic contours drawn arc 
smaller than the statistical contours. 

We will use the functional form filling approach in the 
next Section to place requirements on weak lensing system- 
atics. For general conclusions please skip to Section 6. 



5 AN APPLICATION TO COSMIC SHEAR 
SYSTEMATICS 

We will now use the functional form filling approach to ad- 
dress shape measurement systematics in cosmic shear (due 
to methods Heymans et al, 2006, Massey et al., 2007; or 
PSF inaccuracies Paulin-Henriksson et al., 2008; Hoekstra, 



2004). This Section represents an extension of the work of 
Amara & Refregier (2007b) where certain particular func- 
tional forms were investigated. Here we extend the analy- 
sis to include all functional behaviour to fully address the 
impact of the systematic. For a thorough exposition of cos- 
mic shear we encourage the reader to refer to these exten- 
sive and recent reviews and websites (Munshi et al., 2006; 
Bartelmann & Schneider, 2001; Wittman, 2002; Refregier, 
2003; http://www.gravitationallensing.net). The source 
code related to the work in this Section is released through 
http://www.icosmo.org, see Appendix F for details. 

5.1 Background 

Cosmic shear tomography uses both the redshift of a galaxy 
and the gravitational lensing distortion, shear, to constrain 
cosmological parameters. The observable in this case is the 
lensing power spectrum as a function as redshift and scale 
Ce(z). Since we have redshift information the galaxy popu- 
lation is split into redshift bins where each bin has its own 
auto-correlation function and the cross-correlations between 
bins are also be taken. 

Throughout this Section we will use a fiducial cosmol- 
ogy of Q m = 0.3, Q.DE = 0.7, n B = 0.0445, h = 0.7, w = 
—0.95, w a — 0.0, as = 0.9, n s — 1.0 where we parameterised 
the dark energy equation of state using w(z) — wo+w a (l—a) 
(Chevallier & Polarski, 2001; Linder, 2003) and included the 
spectral index n s , we consider non-flat models throughout 
where Qj, = 1 — S7 m — Qde- We assume a weak lensing sur- 
vey (similar to the DUNE/Euclid proposal, Refregier et al., 
2008a) which has an Area= 20, 000 square degrees, a median 
redshift of z m = 0.8 (using the n(z) given in Amara & Re- 
fregier, 2007b) with a surface number density of 40 galaxies 
per square arcminute. We also assume a photometric red- 
shift error of <r z (z) = 0.03(1 + 2) and split the redshift range 
into 10 tomographic bins. 

The observed lensing power spectrum can be written as 
a sum of signal, systematic and noise terms 

so that an estimator of the observed lensing power spectrum 
can be written by subtracting the shot noise term C" olso 

<5p s = C\ cns + C s / S (33) 

where C( enB is the underlying true lensing power spectrum, 
dependant on cosmology, which we want to measure and 
Cf B is some residual systematic. The error on this estimator 
can be written as 

A Ce = [c\- + cr + crn m 

note that this is the error on the observed signal, not the 
observed signal itself, so that as the systematic increases the 
error on the observed Ce increases. The Fisher matrix and 
bias are then defined in the usual way (equations 5 and 12) 
where oc = AC(. Table 4 shows the expected marginal er- 
rors on the cosmological parameters using the fiducial survey 
calculated using the Fisher matrix formalism, note that no 
prior has been added to these results, they are for lensing 
alone. 

The additive systematic which we consider here is a 
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Parameter 


Marginal Error 




0.006 


^DE 


0.036 




U.U10 


h 


0.086 


w 


0.046 


W a 


0.152 




0.009 


n s 


0.019 



Table 4. The marginal error on each cosmological parameter for 
the fiducial weak lensing survey with no residual systematic. Note 
that no priors have been added on any parameter. 



special kind that has the same shape as the lensing power 
spectrum but where it is multiplied by some unknown sys- 
tematic function (sometimes referred to as a multiplicative 
systematic) Ai (we use the notation of Amara & Refregier, 
2008b) 

Cf s ' y = AiCf (35) 

where ij means the correlation between redshift bins i and 
j, an auto-correlation is where i — j. The particular form 
of multiplicative bias we consider is that which causes the 
true shear, as function of angle and redshift 7 lcns (#, z), to 
be incorrectly determined such that the residual systematic 
shear 'y sys (9, z) is related to the true shear by some function 
m(z) 

^ s (6,z) = m(zW cns (9,z). (36) 

So that the systematic given in equation (35) can be written 
C sys ' lJ = [m(zi) + m(zj)]C l e : ' . This expression from Amara & 
Refregier (2007b) actually assigns an m(z) to each redshift 
slice where each m(z) is the bin weighted average. A more 
accurate approach is to include m(z) in the integrand of the 
lensing kernel (£+/- m Appendix A of Amara & Refregier, 
2007b). We compared results when the bin- weighted aver- 
age was used against the correct integral method and found 
exact agreement, this is because the bins are narrow and the 
functional variation on sub-bin width scales is small. 

We note that Huterer et al. (2006) have used Cheby- 
shev polynomials in weak lensing systematic analysis, al- 
though they marginalise over a systematic parameterised 
using Chebyshev polynomials and then investigate the bi- 
ases introduced if the most likely value of the estimated 
coefficients were incorrect. Here we are not measuring the 
Chebyshev coefficients but rather using this basis set to draw 
every possible function and determine the bias introduced 
by the function itself not some misestimation of any partic- 
ular coefficient. The approach outlined in this article also 
has a resonance with the work of Bernstien (2008) in which 
flexible basis sets are used to address systematic quantities, 
although this is done within the self-calibration (marginalis- 
ing) framework, and their general thesis is that of tuning the 
models that are marginalised over. Bernstien (2008) notes 
that an approach such as the one outlined in this article, 
using the formalism of Amara & Refregier (2007b), would 
be desirable in certain situations. 



N F =10 4 , 0rder=35 

t-H 
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Figure 10. This plot shows the bounded area defined using equa- 
tion (37) with mo = 1 X 10 -3 and f3 = 1.0, shown by the solid 
black lines. This area is then pixelated and the functional be- 
haviour that can occur at each pixel is measured. The colours 
correspond to the percentage of functional behaviour that has 
been experienced for each pixel, this is described in detail in Ap- 
pendix B. 

5.2 Functional Form Filling for m(z) 

We will now place constraints on the function m(z) such 
that that the measurement of the cosmological parameters 
are robust. We first define some boundary (bound function) 
on m(z) such that any systematic function must lie within 
the bounded region. We parameterise the hard boundary 
(Section 4.1) using 

\m(z)\ = m (l + zf , (37) 

we want to know what values of mo and f3 are sufficient to 
ensure that the bias on cosmological parameters b(9i) are 
less than their statistical error a(8i); \b(9i)/a(6i)\ < 1. We 
stress that this is a parameterisation of the boundary within 
which a function must lie, not the function itself. 

To fill the bounded area with every functional form, 
complete down to the scale of Az = 0.2, we use Cheby- 
shev polynomials with a maximum order of -/V = 35 and 
investigate Nf = 10 4 realisations of the basis set, using the 
formalism outlined in Section 4. We do not expect any of 
the cosmological parameters to introduce features into the 
lensing power spectrum on scales Az < 0.2 so this should 
be sufficient. Fig. 10 uses the diagnostic measure described 
in Appendix B to show that on the scale of Az = 0.2 every 
possible functional behaviour has been experienced at every 
point within the bounded region, for this example we use 
m = 1 x 10" 3 and (3 = 1.0. 

To begin we will first consider m(z) with mo = 1 x 10~ 3 
and (5 = 1.0. Using the functional form filling technique we 



© 0000 RAS, MNRAS 000, 000-000 



14 T. D. Kitching et al. 



Parameter 


bias/marginal error 




-0.837 




0.813 


o 


— U.4ol 


h 


-0.483 


Wo 


1.253 


U>a 


-0.999 


T8 


0.855 


n s 


0.868 



Table 5. The ratio of bias to marginal error for each cosmolog- 
ical parameter for the fiducial weak lensing survey and a mul- 
tiplicative systematic of the form given in equation (37) with 
mo = 1 X 10 — 3 and /3 = 1. A negative value implies that the 
bias is negative. 
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Figure 11. This plot shows the bounded area defined using equa- 
tion (37) with mo = 1 X 10~ 3 and j3 = 1.0, shown by the solid 
black lines. Plotted within this area are the functions, out of all 
the possible functions that could be drawn within the bounded re- 
gion, that cause the largest bias on each cosmological parameter, 
denoted by the panel title. 



have identified which systematic functions cause the largest 
bias for each cosmological parameter. Table 5 shows the 
maximum bias to marginal error ratio for each cosmolog- 
ical parameter due to this particular m(z). 

We have taken into account all degeneracies between 
parameters in this exercise - in the terminology of equa- 
tion (12) the inverse Fisher matrix is marginalised over all 
parameters and the Bj takes into account degeneracies be- 
tween the cosmology and systematic functions. 

Fig. 11 shows the worst functions, that cause the largest 
bias, for each cosmological parameter. Since there are degen- 
eracies between all parameters the worst function is practi- 
cally the same for all parameters. 

We find that the ratio of the maximum bias to 
the marginal error is < 1 for most parameters, except 
b(wo)/cr(wo) = 1.25 which is is not acceptable: the bias 



taa. h 




Figure 12. The ratio of maximum bias, found using functional 
form filling, to marginal error as a function of mo and (3 for wo- 
The gray scale represents the bias to error ratio with a key given 
on the side of each panel. The solid lines show the b(fii)/a{Bi) = 1 
contours for each parameter. 



is larger than the statistical error. These results are in 
rough agreement with Amara & Refregier (2007b) where 
it was concluded, with a restricted functional parameter- 
isation of m(z), and the assumption of a flat Universe, 
that mo = 1 x 10~ 3 and (3 — 1 yielded biases that were 
b(8i)/a(6i) < 1.0 but were somewhat on the limit of what 
is acceptable. Furthermore it was concluded that functions 
which have a variation (cross from positive to negative) 
about z > 1 have the largest effect, we also find that the 
functions that cause the largest bias in all the parameters 
varies in the region of z ~ 1. This is because it is at z ~ 1 
that dark energy begins to dominate and so the parameters 
wo and w a , and through degeneracies the other parameters, 
are affected by systematic functions that introduce variation 
at this redshift. 

Fig. 12 shows the ratio of bias to marginal error as a 
function of mo and (3 for wo ■ It can be seen that the redshift 
scaling (3 has the expected effect on the maximum system- 
atic bias: as /3 increases the maximum bias for a given mo 
also increases as the systematic bounded area expands. The 
solid lines in Fig. 12 show the b(0i)/a(6i) = 1 contours. For 
mo = 1 x 10 -3 we need a (3 < 0.70 for the bias on vjq to 
be acceptable. If the redshift scaling is eliminated (3 = 
then the requirement on the absolute magnitude of m(z) is 
relaxed to m <2 X 10~ 3 . 

Fig. 13 shows the affect of survey area on the re- 
quirement of mo and (3. The lines in this figure show the 
\b(wo) /cr (wo)\ = 1.0 contours for varying survey area. We 
have picked too as an example since this parameter provides 
the most stringent constraints on the shape measurement 
parameters (see Table 5 and Fig. 12). As the survey area 
increases and the marginal error on wo decreases the re- 
quirement on shape measurement accuracy becomes more 
stringent. We have fitted a simple scaling relation to these 
contours so that for statistical errors to be reliable the fol- 
lowing relation holds 

2.4 



0.17 



mo 



1 x 10~ 3 



/ Area \ 15 
V 20000 J 



10 p < 1. 



(38) 



The requirements on mo and (3 for the largest survey 
considered are at the limit of currently available shape mea- 
surement techniques, that yield on average m ~ 2 x 10 -3 
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Figure 13. The solid lines show the \b(wo) /o-(wq)\ = 1.0 contours 
in the (mo, 0) parameter space for varying survey area. Black 
(solid) shows the contour for the fiducial survey with Area= 20000 
square degrees, red (light gray, dashed) shows the contour for a 
survey exactly the same as the fiducial survey expect that Area= 
2000 square degrees, and green (lightest gray, dot-dashed) for a 
survey with Area= 200 square degrees. For the uifj constraint to 
be robust to m(z) systematics the values of mo and (3 must lie 
leftward of the contours, see Fig. 12. 

at best (Miller et al, 2008; Kitching et al., 2008b). A large 
redshift variation of m(z) results in large biases (Fig. 11) so 
a shape measurement method that is unaffected by the mag- 
nitude/size of the galaxy population should yield robust cos- 
mological constraints. For example Kitching et al. (2008b) 
have shown that the lensf it method has a characteristic m 
that has a small magnitude/size dependence. In this case, 
<C 1.0, and a more relaxed constraint on mo ~2x 10~ 3 - 
4 x 10 -3 is required which is well within reach of these most 
recent developments in shape measurement. These results 
are also in agreement with Semboloni et al. (2008) where 
they find that for a very restrictive class of m(z) functions, 
of the form m(z) — az + b, but a more complex likelihood 
description, that some parameter combinations of a and b 
can give rise to biases that are less that the cosmological 
errors. 

For the smaller surveys considered, that are well 
matched to currently available or upcoming experiments 
(e.g. CFHTLS, van Waerbeke et al., 2001; Pan-Starrs, 
Kaiser, 2002), the shape measurement techniques currently 
available have biases m (Heymans et al., 2006; Massey et 
al., 2007; Kitching et al., 2008b) that are well within the 
required level of accuracy. The caveat to these conclusions 
is that here we have not discussed the size or galaxy-type 
dependence of any bias and we have not considered any cal- 
ibration offset in the measured shear. 



6 CONCLUSION 

In this article we have presented a method that allows one 
to move beyond the tendency to treat a systematic effect 
as a statistical one. When a systematic is treated as a sta- 
tistical signal extra parameters are introduced to describe 
the effect, and then these extra nuisance parameters are 
marginalised over jointly with cosmological parameters. We 
have shown that even in a very simple toy model such an 



approach is risky at best. The results being highly depen- 
dent on the choice of parameterisation, and in the limit of a 
large number of parameters any statistical signal on cosmo- 
logical parameters can be completely diluted. One could use 
marginalisation if there exists a compelling physical theory 
for the systematic, however if the functional form is merely 
phenomenological or if it contains a truncated expansion 
then one should use caution. Even in the case that confi- 
dence is high with respect to the assumed functional form 
any residual must be investigated correctly. 

As an alternative we advocate treating a systematic sig- 
nal as such, an unknown contaminant in the data. We ad- 
dress the situation where we have at least some extra in- 
formation on the systematic, either from some theory or 
simulation that may provide a hard boundary within which 
a systematic must lie, or from some external data set. We 
leave the case of 'entangled systematics' (where there is no 
extra information and where the systematic depends on the 
cosmological parameters) to be investigated in Amara et al. 
(in prep). Throughout we have introduced each concept us- 
ing a simple toy model. 

Since the systematic is treated as a genuine unknown 
we must address every possible functional form that is al- 
lowed, by either the hard boundary or the extra systematic 
data. To do this we use complete basis sets and randomly 
sample from the space of coefficients until all functional be- 
haviour has been experienced at every point with the hard 
boundary, or within a few a of the data. We have shown that 
such "functional form filling" can be achieved by this tech- 
nique using either Chebyshev, Fourier or tophat basis sets. 
By treating each function as a possible systematic, a bias 
in the maximum likelihood value is introduced whilst the 
marginal error stays the same. Throughout we have shown 
that as long as functional filling is achieved all conclusions 
on the magnitude of a systematic effect are independent of 
the choice of complete basis set - though binning is highly 
inefficient in terms of computational time, and could be la- 
belled as an unphysical basis set. 

For the case of a hard boundary every function is given 
an equal probability and as such a maximum bias exists. For 
the case of extra systematic data we show that a probability 
can be assigned to each function, and bias, allowing for the 
production of robust systematic probability contours. 

We have made a first application of hard boundary func- 
tional form filling by addressing multiplicative systematics 
in cosmic shear tomography; we have left an application of 
the data bound for future work. We address a lensing sys- 
tematic, that can result from shape measurement or PSF 
reconstruction inaccuracies, that has the same overall shape 
as the lensing power spectrum but is multiplied by some 
extra function. This is commonly represented using a multi- 
plicative function m(z) (Heymans et al, 2006; Massey et al., 
2007) where 7 sys = m(z)j tIuc . For a DUNE/Euclid type sur- 
vey, we find that in order for the systematic on wo to yield 
a bias smaller than the marginal error the overall magni- 
tude of 771(2) should be mo < 8 x 10 -4 for a linear scaling in 
(1+z), but as the magnitude of the redshift scaling is relaxed 
then the requirement on mo increases to mo < 2 x 10 -3 . The 
most recent shape measurement methods have m~2x 10~ 3 
(e.g. using the lensf it method Miller et al., 2007; Kitch- 
ing et al., 2008b) and have a small scaling as a function of 
magnitude/size. The results shown here then, coupled with 
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some currently available shape measurement techniques im- 
ply that constraints on cosmological parameters from to- 
mographic weak lensing surveys should be robust to shape 
measurement systematics (with the caveat that PSF cali- 
bration and galaxy size dependence of the bias have been 
neglected here). 

The techniques introduced here should have a wide ap- 
plication in cosmological parameter estimation: in any place 
in which there is some signal with some extra information on 
a systematic. For example baryon acoustic oscillations and 
galaxy bias, CMB and foregrounds, galaxy clusters and mass 
selection, supernovae and light rise times, weak lensing and 
photometric redshift uncertainties. A sophistication of these 
techniques could assign different weights/prior probabilities 
to particular functional forms, that are known a priori to 
have a large/small effect on cosmological parameter deter- 
mination, such weights could come from either theoretical 
or instrumental constraints. 

Cosmology is entering into a phase in which the sta- 
tistical accuracy on parameters will be orders of magnitude 
smaller than anything achieved thus far. But we must take 
care of systematics in a rigorous way to be sure that our 
statistical constraints are valid. There exists the pervading 
worry that the functional forms used to parameterise sys- 
tematics are not representative of the true underlying nature 
of the systematic and that something may have been missed. 
In such a scenario one should always add the warning that 
"cosmological constraints are subject to the assumption of 
the systematic form" . 

Here we have presented a way to address systematics 
in a way that requires no external assumptions, allowing for 
robust and rigorous statements on systematics to be made. 
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written for the joint log-likelihood of the cosmological pa- 
rameters and the systematic 

2£(0, S)=Y. a c 2 {C-T- Sf + J2 °d 2 (D - S) 2 . (45) 



APPENDIX A : MARGINALISATION VS. BIAS 

Here we will show how marginalisation is mathematically 
different to the bias formalism used in the main article. 

When marginalising the log-likelihood of some cosmo- 
logical parameters can be written, for a signal C some 
theory T and a systematic S as 



(39) 



When marginalising over the systematic the effect is pa- 
rameterised by some extra parameters a. If extra data D is 
provided on the systematic then this adds a prior such that 
the total log-likelihood can be written 

2C(0, a)=J2 a c 2 (C-T- S) 2 + £ a D 2 (D - S) 2 . (40) 

X X 

The Fisher matrix for this is created by taking the deriva- 
tives of the log-likelihood with respect to cosmological pa- 
rameters 6 and systematic parameters a so that the first 
and second terms in equation (40) give the following Fisher 
matrices 



F° 



F ea 

?aa ] + 

c 
7 ea 





paa 



7 ae 



(41) 



The marginal error on the cosmological parameters is found 
by inverting equation (41). The inverse of the upper left- 
hand segment of the Fisher matrix in equation (41) can be 
found by using the Schur complement of the block matrix 
and then expanding this using the Woodbury matrix iden- 
tity which gives 

(-f **r-icft) 1 = MSE marg 

= A' 1 + A~ 1 B(E - B T A' 1 B)' 1 B T A' 1 

(42) 

where A = F 6e , B = F ea and E = F§ a + Fg a . The 
marginal error in the case of no systematic (A^ 1 ) has been 
increased by an extra factor that depends on the degen- 
eracies between the systematic parameters and the cosmo- 
logical ones (B) and on the information available on the 
systematic parameters themselves (E). This is also equal to 
the mean square error of the cosmological parameters in this 
case, since no bias is introduced. 

For the bias functional form filling technique the log- 
likelihood of the cosmological parameters can again be writ- 
ten as 



2C(0)=J2 a c 2 (C- T -S) 2 



(43) 



In the data bound case we can assign a weight to each sys- 
tematic function S (not extra parameters) 

2£(S)=J2 a o 2 (D-S) 2 (44) 

X 

so that a very similar expression to equation (40) can be 



We have not necessarily parameterised S with any parame- 
ters we wish to measure, S can simply be a function that has 
been arbitrarily drawn through the systematic data. How- 
ever if the systematic is parameterised (as we do in the func- 
tional form filling approach using complete basis sets) then 
the signal data C is not used to determine the values of the 
extra parameters - the number of degrees of freedom in the 
fit has been reduced with respect to the marginalising case. 
The systematic data D (or a hard boundary) is used to de- 
termine the probability of the systematic. Hence the Fisher 
matrix can be written, in this case, like 







1 




(I 


F aa J 







(46) 



The inverse of the upper left-hand segment is now simply 
(F 96 )^ 1 = A~ l . Using the result from Appendix D (for the 
hard boundary case) a bias is introduced with a maximum 
value of 



max|6i| = (const) (.A 



90-. 



= (const) (A 



(47) 



Fj = 



We have again condensed the notation so that 

~ 2 f^. Note that F ^ B ^ E. Hence the total mean 



J2°c z w-- Note that F + B + 
square error on cosmological parameters for the bias case 
comes from the unaffected marginal error of the cosmologi- 
cal parameter and the bias 

MSE bias = A' 1 + bias 2 = A' 1 + (const) 2 (A' 1 ) 2 F 4 (48) 

a similar expression can be written for the data bound case, 
see Appendix D equation (85) . We have compressed the sub- 
script notation from equation (47) in equation (48). 

These equations show how marginalisation and the bias 
formalism relate in a general sense, and that they are math- 
ematically different objects i.e. when finding the maximum 
bias using the functional form filling approach the MSE be- 
tween marginalisation and functional form filling can be dif- 
ferent. This difference arises because the space of functions 
probed by marginalisation is restricted. In the next section 
we will show that if the functional space assessed is the same 
then the MSE should be equal. 

Information Content 

If the parameterisation, and number of (truncated) param- 
eters, are the same then the MSE resulting from the bias or 
marginalising is the same. We show this here using a simple 
illustrative example. 

We approximate the nuisance correlation C sys as a first- 
order expansion of a set of parameters 



C By » 



CJ 5 ™ + aV a C° y 



(49) 



This can be though of as a restricted set of functions or a 
special case of C sys . The nuisance parameters a will be cor- 
related with the cosmological parameters 8 and we can form 
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the extended vector and Fisher matrix as shown in equation 
(41), as in equation (42) we can write the marginalised co- 
variance matrix of like 

{60 t )c = [F<S,4,]gg 

= F g -/ + F^FeaiFaa - Fg a Fg~ g 1 Fg a y 1 Foa FJffi ) 

and there is no bias of the measured parameters. The sub- 
script 'c' on the covariance indicates that we are including 
the covariance between and a. 

If instead we do not account for the covariance between 
and a in the data, the measured error in the parameters 
is 

(90'} = [Fgg]- 1 (51) 

but we have induced a bias in the measurement of 8, 

„ ^/""signal 

A0 = -F^£^ 2 ^-C sy >). (52) 

where a 2 D is the variance of the data. Using our linear ex- 
pansion of C sys (a) we find 

A0 = -Fg^Fgaa. (53) 

From this can see the bias effect explicitly comes through the 
correlation of the nuisance parameters with the cosmological 
parameters, via Fg a . 

We can simplify things further by considering a single 
parameter and a single nuisance parameter a. The condi- 
tional error on the parameter can be simplified to 

(00) c = {60){1 - r 2 )" 1 (54) 

where we have introduced the correlation coefficient defined 

as, 



r 2 — Fgg 1 Fg a F a a ■ 



(55) 



We can see this by looking at the inverse of the 2x2 matrix 



F&& — 



J / F aa 

- F 6a V ~ F "< 



-Fg a 
g FggF a , 



(56) 



FggF a 

The marginalized error on 9 is 

(99) c = f aa =(99)(l-r 2 )- 1 - (57) 

reeraa r 6a 
Similarly the marginalized error on a is given by 



(aa) c = (aa)(l — r)~ 



(58) 



Note that the marginalized errors on 9 and a go to infin- 
ity when the correlation coefficient is unity, r = 1. This is 
because 9 and a are completely degenerate. 

From the definition of the correlation coefficient we can 
write the covariance between 9 and a as 



(9a)l = (69) c (aa) c r 2 



(59) 



which agrees with the definition of r given in equation (55). 

Note that the results for (00) c , (aa) c and (0a) 2 are not 
just the errors and covariance when we marginalize. They 
are also the errors and covariances when the parameter 9 is 
correlated with a. The only way to get back the uncorrelated 
errors is if 9 and a are uncorrelated, so r — 0, or if we 
know a from some other measurement, in which case they 
effectively decorrelate. Even if we assume some value, or 



range of values, for a, there is still a correlation between 
possible values of a and 9. 

Now consider the bias on 9, AO, when we ignore the 
nuisance effect (or assume some fixed form) . In our simple 
model the bias can be written 



AO 



1/2 



(99) 



1/2 



<a«> 1/2 (aa)V 2 



(60) 



Now the error on the bias, including all the correlations be- 
tween parameters, is given by 

( (aa),r 2 = (09), r 2 = (00) 

" (61) 

One might have mistakenly thought that the covariance be- 
tween 9 and a was 



(A9A9) C = (aa) c r 2 = (09) c r 2 

\(aa) c 



(A0AO) C = 



(aa) 



(aa)r = 



(62) 



but in doing so we have ignored the real correlation between 
and a that exists. In the case that of fully correlated pa- 
rameters, when r — 1, the error on A9 is only the uncorre- 
lated error on 9, whereas the real uncertainty is infinite. 

Finally the correlation between the bias and the cosmo- 
logical parameter is 



(0A0) C 



1/2 



(aa) 



1/2 



r(0a) c = 



1/2 



\V2/ 



(aa) 



1/2 



(aaYJ 2 r = (00) c r z . 

(63) 



Now we want to compare the effect of the bias with 
marginalization on the error of 0. In the case of marginal- 
ization this is given by the correlated covariance 



(09) c = (00)(l-r 2 )- 1 . 



(64) 



In the case of bias, we have to add the uncorrelated error, 
(99), with the error in the bias value, (A0A0) c , which does 
include the correlation between cosmological and nuisance 
parameters, to form the MSE; 

MSE = (00) + (A0A0) c 

= (00) + (00) c r 2 = (00) + (00) 



(^(l-r 2 )- 1 



(65) 



Hence the MSE is the same as the marginalized error in this 
simple case, and there is no loss or gain of information. Our 
result here is for a 2-parameter case, but is easily extended 
to multiple parameters. 

It is worth thinking about the assumptions that lead to 
this result. For instance we have assumed in our analysis that 
the nuisance bias is completely specified by the parameter a. 

This shows that if one fully understands the nuisance 
effect, there is no difference between marginalisation and 
form-filling. Given marginalisation can be done quickly, and 
in some cases analytically, we would advocate marginalisa- 
tion in this case. 

However, if the form of the nuisance function is not 
known, or is wrong, marginalization will tend to underesti- 
mate the true systematic error because of the limited func- 
tional space explored by the set of functions assessed by the 
truncated basis set. 

Given that in general the number of parameters that 
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can be marginalised over is limited by the number of data 
points if the systematic parameterisation is unknown then 
the marginalising will necessarily not be able to assess the 
impact of some functions. 



Constraining nuisance effects with external data 

The effect of additional, external data to constrain the nui- 
sance effects will add an extra Fisher matrix to F aa ', 



F — > F' — F 4- F c 

1 aa f 1 aa — ± act T J a 

If we define a new parameter 

P = [-faa] -^aa 



'7 C 



(aa) e 



(66) 



(67) 



as the ratio of the conditional error on a from the original 
dataset to the expected measured accuracy on a from the 
external data. If the error on a from the external data is 
small, we can expect (3 to be large, while of the external 
error is large, j3 will be small. This results in the correlation 
coefficient becoming 



(68) 



l + P 

When the error on the nuisance parameters is highly con- 
strained by the new data set, > 1 and the correlation 
coefficient decreases. The effect of new data in the nuisance 
functions is to de-correlate the nuisance parameter from the 
data. Propagating this through we find 



(aa)' c = (aa) ( ^—±—^ j 



(69) 



so that the improved accuracy from the external dataset 
feeds through. The marginalised error on the cosmological 
parameter is now 



1 - r 2 + (5 



The error on the bias is now 

(A0A0) = (eeyy 

Finally, the MSE is 

MSE' = {99)' + {A9A9)' C 



(l-r 2 + /3)' 



(70) 



(71) 



(1 



r 2 + j3) 
■p 



1 



+ P 



(72) 



Hence adding an external dataset to constrain a does not 
change the conclusions that, if the systematic effect can be 
modelled by a known parameterization, the MSE is the same 
as the marginalized error on 9. 

Summary 

We summarise this Appendix by stating the its main con- 
clusions 

• If the functional form of the systematic is known then 
the degradation in the MSE as a result of marginalising or 
functional form filling (bias) is the same. 
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Figure 14. For each pixel defined within a bounded region we 
label the neighbouring pixels. The function shown would assign 
the functional dependency (enter 2, exit 6) to the central pixel. 



• If the functional form of the systematic is unknown then 
the MSE from marginalising will tend to underestimate the 
true systematic error and is in general not equal to the MSE 
from functional form filling. 

• Given that the marginalising necessarily truncates the 
basis set there are always some functions that marginalising 
cannot assess. 



APPENDIX B: TESTING FUNCTIONAL 
FILLING 

In this Appendix we will present a numerical analysis that 
will show that if the order and the number of function evalu- 
ations is sufficient then a bounded area can be exhaustively 
filled with all possible functional forms complete to some 
scale. We want a non-parametric minimum-assumption ap- 
proach for determining whether we have fully sampled the 
allowed function-space. 

We can pick a certain scale upon which we can inves- 
tigate whether all possible functional forms are evaluated. 
Having chosen a scale a bounded area can now be pixe- 
lated at that scale. In Fig. 14 we consider a particular pixel 
and its neighbouring pixels. There are 22 non-degenerate 
ways in which a function can pass through the pixel in ques- 
tion and two surrounding pixels (entering via one pixel and 
exiting via another) for example (enter 1, exit 2), (enter 1, 
exit 3), (enter 1, exit 4), (enter 8, exit 6). We call each 
combination of entrance and exit a 'functional dependency' 
or a 'functional behaviour'. We exclude functions that are 
multiple- valued e.g (enter 8, exit 7) which would require 
multiple values of y(x) for a given x. 

Some combinations are more likely than others, for ex- 
ample consider a box centered on (0,0). The function y = x 
will enter through 1 and exit via 5, corner functions of this 
type are unlikely but not impossible; a more robust and fair 
measure may concatenate boxes 1 and 2 (and 3 and 4) for 
example. 

We now look at each and every pixel within the 
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Percentage of Functional Forms Filled 

Figure 15. For the Chebyshev basis set. For each pixel de- 
fined within the simple bounded region we sum up all functional 
dependencies experienced from the full set of functions evaluated. 



20 40 60 80 100 

Percentage of Functional Forms Filled 

Figure 16. For the Fourier basis set. For each pixel defined 
within the simple bounded region we sum up all functional de- 
pendencies experienced from the full set of functions evaluated. 



bounded region and for each function drawn perform a check 
of which entrant and exit functional dependencies are ex- 
plored - checking off each functional dependency as it is 
sampled. We then assign for each pixel the percentage of 
functional dependencies that have been experienced as a re- 
sult of having drawn the full set of functions. We perform 
this test for the hard boundary toy model presented in Sec- 
tion 3.3 for each of the basis sets considered in Section 4, 
Chebyshev, Fourier and tophat functions. We increase the 
maximum order N in the expansion in equation (22) and the 
number of random realisations Nf of the coefficient space 
{ao,...,ajv} (and {&o,-..,&iv} for the Fourier basis set) . 

Figs 15 to 17 show the toy model hard boundary that 
has been pixelated on the scale of Ax = 0.5. It can be seen 
that within the bounded area every pixel does experience all 
functional behaviour if the order and number of functions 
drawn is sufficiently large. As expected when the functional 
order increases, and as the number of function evaluations 
increases, the number of pixels that experience all functional 
dependencies increases. We find that a maximum order of 
N >35 with Nf > 10 4 realisations is sufficient for each point 



in the bounded area to experience every functional form for 
this simple example. 

For the tophat basis set (binning) we find that the max- 
imum order (number of bins) needs to be much larger than 
for either Chebyshev or Fourier basis sets. This can be un- 
derstood since if the bin width (order) is larger than the 
resolution of the pixels then it is impossible for any given 
pixel to experience particular functional dependencies (such 
as [enter 1, exit 2]). So even though tophat functions can be 
used for functional form filling there is some computational 
expense in this choice. 

Fig. 18 shows some diagnostic plots representing the 
completeness of the functional space sampled. We show the 
percent of pixels that have 100% of the possible functional 
behaviours sampled, and the percentage of functional be- 
haviours sampled by the average pixel. In the top panels we 
vary the scale that is investigated, and it can be seen that 
for the Chebyshev and Fourier basis sets ~ 100% of pixels 
have experienced every functional behaviour down to a scale 
of Aa; ~ 0.5. The tophat basis set performs much worse with 
a mean filling of ~ 10%, and only ~ 50% of all pixels expe- 
riencing every possible behaviour. The general trend is that 
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Figure 17. For the tophat basis set. For each pixel defined 
within the simple bounded region we sum up all functional de- 
pendencies experienced from the full set of functions evaluated. 
For the tophat basis order corresponds to 'number of bins'. 



as the scale drops below the average oscillation length of the 
most highly varying functions the percentage of 'good' pix- 
els (experiencing all functional behaviour) sharply declines. 
The tophat basis set fails at Aa; ~ 0.3 since, with N = 35 
this is approximately the 'bin width' of the functions. In a 
cosmological application one would ensure that the form fill- 
ing was complete down to the scale upon which a particular 
parameter affects the signal. 

In the middle panels of Fig. 18 we vary the maximum 
order of the basis set expansion whilst keeping both the 
scale and the number of random realisations fixed. It can 
again be seen that for a maximum order > 35 the Chebyshev 
and Fourier basis sets achieve 100% functional filling. Again 
the tophat basis set (binning) fails to achieve any complete 
functional filling, and only begins to fill in some pixels when 
the number of bins becomes more than the number of pixels. 

The lower panel of Fig. 18 shows how the filling effi- 
ciency varies with the number of random realisations of the 
basis set. For the simple example given the Chebyshev and 
Fourier basis sets completely fill the bounded area, down 
to the scale of Ax — 0.5 in ~ 10 3 realisations whereas the 
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Figure 18. The right panels show the percentage of pixels in the 
bounded area that have experienced every type of functional be- 
haviour (100% filling). The left hand panels show the percentage 
of all functional behaviours experienced per pixel on average. The 
upper panels shows these diagnostic measures of functional form 
filling as a function of scale, with the order and number of reali- 
sations fixed at N = 35 and Np = 10 3 respectively. The middle 
panels as a function of maximum order, with the scale and num- 
ber of realisations fixed at Ax = 0.5 and Np = 10 3 respectively. 
And the lower panels as a function of number of random reali- 
sations, with the scale and order fixed at Ax = 0.5 and N = 35 
respectively. In all panels blue (darkest gray, lower lines) is for the 
tophat (binning) basis set, red (lighter gray) for the Chebyshev 
basis set and green (lightest gray) for the Fourier basis set. 

tophat basis set requires many orders of magnitude more 
realisations. 

One could imagine futher metrics that could be used 
to gauge the completeness of the functional space that has 
been sampled, in this first exposition of the methodology we 
have shown a simple way to gauge this effect. 

One concern is that this metric may favour the tophat 
basis since we are using a pixelated measure of scale. How- 
ever we find that even with this advantage the tophat ba- 
sis set achieves complete filling only at the expense of a 
prohibitive amount of computational time compared to the 
Chebyshev and Fourier basis sets, as the order and the num- 
ber of realisations would need to be much larger to compen- 
sate for the pathological choice of basis functions. 

We note that a step-function (e.g. function / over the 
interval x £ [0, 1] which is equal to +1 for < x < 0.5 and —1 
for 0.5 < x < 1) may be better approximated by a low order 
tophat function and may take a very high order Chebyshev 
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or Fourier expansion. The feature highlighted here - a step 
function - has a feature (the step) which has a very small 
scale variation. In fact the scale of a step is actually zero. 
The issue of wether functional form filling is complete down 
to some scale is shown in Fig. 18, and as the maximum 
order is increased and more highly oscillatory functions are 
included the minimum complete scale will decrease. This 
example highlights that fact that the tophat basis set can 
probe these very small scales to some degree whilst missing 
larger scale features, this can be seen in Fig. 17. 

For all the calculations in Section 4 we pessimistically 
use a maximum order of N = 35 and use Nf — 10 4 realisa- 
tions for all basis sets. 



APPENDIX C : COMPUTATIONAL TIME 

One may be concerned at the amount of computational time 
that the functional form filling we are advocating may take. 
This is a valid convern since many realisations of the sys- 
tematic data could be required and one must re-perform the 
cosmological parameter fitting for each and every possible 
function in order to evaluate the potential bias. 

If, for a simple grid search in parameter space, the time 
taken to analyse a single point in parameter space is r 
then for N cosmological parameters the total time scales 
as T oc A N t where A is the number of evaluations per pa- 
rameter. Even for one parameter A > 1 to correctly map 
out a likelihood surface. If one marginalises over M extra 
systematic parameters then the total time must increase to 
Tmarg oc A (n+m ^t. For the functional form filling case the 
total time taken for the calculation is simply the number of 
different function evaluations F multiplied by the number 
of data realisations D and the time to estimate the cosmo- 
logical parameters is Tgf oc FDA n t. The ratio in the com- 
putation time can now be written as T marg /Tfff « A M / FD. 

Let us pessimistically assume that approximately F = 
10 4 function evaluations are needed and D = 10 3 data real- 
isations, and conservatively assume that M = 10 extra sys- 
tematic parameters are required to ensure a systematic effect 
is correctly determined. For A > 100 we find T marg > 10 9 Tfff . 
So on the contrary to such a technique being computation- 
ally expensive, treating systematics in such a fashion may 
in fact be much more efficient than marginalising over many 
nuisance parameters using a traditional grid search. 

For some Monte Carlo parameter searches the compu- 
tational time can be reduced to T = NAt. In this case we 
find that T marg /T fff w M/NFD + 1/FD. The computation 
time is now longer for functional form filling T marg /TfFf <C 1. 
This highlights the paramount importance of finding effi- 
cient form filling functions such as Chebyshev polynomi- 
als; we investigate the efficiency of Chebyshev, Fourier and 
tophat basis sets in Appendix B. 

We take random realisations of the coefficient-space to 
uniformly sample the coefficients' possible values. This is a 
first, brute-force attempt at the problem. Alternative ap- 
proaches could explore a set grid of values in this space, or 
use a more sophisticated random search such as a Monte 
Carlo chain approach - we leave this for future work. 

In this first investigation we have taken a simplistic ap- 
proach and, as shown in Appendix B, we have found that 



a set of simple uniform random realisations is sufficient to 
explore the full space of functions. 

Bayesian vs. Frequentist One may be concerned that 
what we are advocating is a frequenist solution to the sys- 
tematic problem. On contrary what we suggest is explicitly 
Bayesian : the method can take into account any prior in- 
formation on the systematic from either theory or data. The 
specific method used in this article to find the form filling 
functions is analogous to a brute force parameter in maxi- 
mum likelihood rescontruction. 



APPENDIX D : EXTREMAL OF THE BIAS 

Here we will show that given a hard boundary or some data 
there exists a maximum absolute bias. Within the Fisher 
matrix formalism the bias caused by a function can be writ- 
ten (equation 12) as 

a^-*si g nal 

b(0i) = (F-% £ s c 2 C SYS ^qq- (73) 

where sc is the observed signal variance. We can rewrite 
this as 



where D a%l , = , we compress this to 



(74) 



(75) 



where = (F 1 ) flv s a 2 D at „. 

We have the additional constraint that we are consid- 
ering the set of systematic functions {S} that give the same 
weight 



(T a 2 (Ca ys — d a ) 2 = A 2 = constant 



(76) 



with respect to some data vector d a . In this proof and 
throughout the article we assume that there is at least some 
systematic information at every data value in te signal. If 
there was no constraint at the position of some of the signal 
data values then the bias would be unbounded. 

For the hard boundary there is a constraint that 
[C sys {x)\ 2 < A 2 at all x - this is actually for a constant (flat) 
hard boundary for a variable boundary the A — > A(x). 

So to find the maximum bias we need to solve the fol- 
lowing equation 



d 



OCT 



\cTQ^ - a (X>" 2 (cr - d a ) 2 - A^j | = 







(77) 

where A is a Langrange multiplier. Solving this for fixed fi 
and all 7 we find that 



C 7 ^ G ~iQi\ii + d-f. 



(78) 



To determine the value of the Lagrange multiplier we sub- 
stitute this back into the constraint equation (76) to get a 
quadratic in (1/2A) 



0. (79) 
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This has solutions of the form 
1 

2A 



where 



~ 2E Q (^) 2 



but would have a smaller bias than a "good function" that 
went through the same data points because the ao behaviour 

(80) 

is not strongly degenerate with some "bad behaviour" but 
has a simple functional dependency. 

The key conclusion of this Appendix is that out of the 
space of all functions that have the same weight with respect 
to the data there exists a maximum and a minimum bias 
that this space of functions can cause. 



2 J2a (cSQ^m) 

So the function(s) with the maximum bias can be expressed 
as 



C S J B = {R^±T, 



+ d~, 



(82) 



Note that there exists two functions that represent the max- 
imum and minimum bias, which are given by 

max/min(& M ) = [R^oiQ^u. + d a ] Q a;(1 ± [T M cr 2 Q a;M ] Q Q:M . 

(83) 

So for a set of systematic functions that give the same 
weight with respect to some data there should exists a maxi- 
mum and a minimum bias. Equation (76) can be understood 
by looking at Fig. 7. For the data boundary case there will 
be some large and small biases - actually every bias is 'al- 
lowed' but has some probability with repsect to the data. 
Equation (76) effectively takes a horizontal cut across on of 
the scatter plot panels in Fig. 7 at a given weight so that 
equation (83) gives the minimum and maximum biases for 
that weight. 

If the mean of the data is zero (d a ) = and the variance 
of the data is (d 2 a ) = a 2 , then {R^} = and 



(^SQ«;m) 



(84) 



(81) APPENDIX E: DATA WEIGHTING 

We know that given some error bars on data, and a different 
realisation of the experiment, that the actual data points 
would be scattered differently but the statistical spread of 
the data would be the same. If a line fits exactly through 
some data points it is given a likelihood of exactly 1 but we 
know that this probability is spurious since given another 
realisation of the experiment the data would be differently 
scattered and a new function would be assigned a probability 
of 1. This is a problem because P = 1 should be unique. This 
problem occurs because in the original step the probability 
of the data has not been taken into account. 

For the specific Gaussian case only we have used the 
marginalisation over the probability of the data can be done 
analytically. We stress here that in general the systematic 
mean will be non-zero and that only in this very simplified 
Gaussian case (and some other analytic examples) can this 
procedure can ve done anlytically. 

Here we will outline how this 'sample variance effect' - 
that of taking into account the likelihood of the systematic 
data - can be incorporated into the \ 2 weighting scheme 
given in equation (26). 

We can write the probability of a function f[a] as the 
sum over the data vector d t multiplied by some prior prob- 
ability of the data values p(di) 



so that the mean bias is zero and the bias contours (contours 
are drawn for a constant weights, A) can be written as 



(i> M ) = ±(T , m) " 2 Q 2 « ; m- 
For a hard boundary we have 



Ru 







T„ = 



A 



= A' 



which gives solutions for the maximum bias 
max/min(6 M ) = ±A'Q 2 a 



) 2 



(85) 



(86) 



(87) 



this calculation can be compared to equations (84) and (85). 
This confirms that in the hard boundary case the absolute 
value of the bias has a maximum and that there exists two 
mirrored functions (f(x) and —f(x)) which both yield this 
maximum absolute bias. 

In the case when the systematic is constrained by some 
data points C s ys a , we can easily construct a function that 
goes through all the systematic data points but has un- 
bounded "bad effects" at some other point. However such 
a function, that has the same fit to the data but some other 
behaviour between the points is also assessed in terms of its 
impact on the cosmological parameter in question. Looking 
at Figure (7) such a function would have the same weight 



P(f[a]) = i[p(fi\di)p(d i 



where /; = /(a;;; a) are the function values are the variable 
positions Xi at which the data has been taken. 

We now want to marginalise over the probability of each 
data point to obtain p(f[a]) 



p(f[a}) = Jp(f[a}\D)dD 

= n/ p(/ii*)p(*)d*- 



(89) 



Since we assume Gaussian data throughout the probability 
of the data can be written as a Gaussian centered about zero 
and p(di\fi) can be replaced with a \ 2 distribution such that 



p(/N) 



„■ J — c 



(dj- 



72 — 



>D dd. 



(90) 



where Oi is the error on the i th data point and and ao is 
the variance of data which we take to be ao = 0». Note 
that each data point could have a different error bar - we 
have made the assumption that the data at each point is 
Gaussian distributed not that all data points are the same 
(the integral is over di not over i). Evaluating the integral in 
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equation (90) we have (now using log-likelihood for clarity) 



So that in the case of Gaussian distributed data the proba- 
bility of each function (and hence the probability of the bias 
incurred by that function) can be simply evaluated using 
equation (91). Note that a Gaussian with a mean of zero 
is unqiuely defined by its variance hence the data values 
themselves do not appear in the weighting formula given. 



APPENDIX F : ICOSMO MODULE 
DESCRIPTION 

The 'worst bias' calculations presented in Section 5 were 
done using an extension to the open source interac- 
tive cosmology calculator iCosmo (Rcfregier et al., 2008b; 
http://www.icosmo.org). This additional module will be 
included in vl.2 and later. 

The tomographic lensing module mk_bias_cheb takes 
mo and f3 defined in equation (37) and uses the functional 
form filling technique (using the Chebyshev basis set) to cal- 
culate the maximum bias in each cosmological parameter. 
The lensing survey and central cosmology can be arbitrarily 
defined using the common set_f iducial routine described 
in Refregier et al. (2008b). 




(91) 
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